import pandas as pd
import geopandas as gpd

500m Grids

grid_path=r"C:\Users\User\Documents\EHI\Projects\VAR Modelling\Data Mgmt\500 grids\grid_500m_clipped.shp"
grid_layer = gpd.read_file(grid_path)

grid_layer.head()

	OBJECTID	Id	Shape_Leng	Shape_Area	geometry
0	1	0	2000.000000	250000.000000	POLYGON ((18162.733 31755.775, 17662.733 31755...
1	2	0	2000.000000	250000.000000	POLYGON ((23162.733 42755.775, 23662.733 42755...
2	3	0	2000.000000	250000.000000	POLYGON ((37662.733 40255.775, 37162.733 40255...
3	4	0	2961.516312	226945.081982	MULTIPOLYGON (((28662.733 48255.775, 28235.903...
4	5	0	1974.913982	248059.096956	POLYGON ((29162.733 28255.775, 29162.733 27755...

grid_layer=grid_layer.set_crs(epsg="3414")

grid_layer.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 3499 entries, 0 to 3498
Data columns (total 5 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    3499 non-null   int64   
 1   Id          3499 non-null   int64   
 2   Shape_Leng  3499 non-null   float64 
 3   Shape_Area  3499 non-null   float64 
 4   geometry    3499 non-null   geometry
dtypes: float64(2), geometry(1), int64(2)
memory usage: 136.8 KB

DBF with XY

df = pd.read_csv("DEN-Hexagon DBF 2013_2020_complete.csv")

df.head()

	CASE_ID	CLUSTER_ID	Onset_EYear	Onset_EWeek	Clean_Onset	ONSET_DATE	NEA_ONSET_	NOTIFY_DAT	POSTAL	X output	Y output	cell_id	friendly_n	Shape_Leng	Shape_Area
0	DENF-2013-000025	0	2012	51	17/12/2012	17/12/2012	17/12/2012	31/12/2012	650536	18857.18765	37512.04877	1334	HONG KAH NORTH	0.016238	0.000019
1	DENF-2013-000018	20130008	2012	51	21/12/2012	31/12/2012	21/12/2012	31/12/2012	536467	32943.97740	36955.30429	2826	UPPER PAYA LEBAR	0.016238	0.000019
2	DENF-2013-000048	0	2012	51	22/12/2012	30/12/2012	22/12/2012	3/1/2013	550421	32663.67780	36774.03736	2826	UPPER PAYA LEBAR	0.016238	0.000019
3	DENF-2013-000109	0	2012	51	22/12/2012	3/1/2013	22/12/2012	4/1/2013	550422	32750.05966	36774.55328	2826	UPPER PAYA LEBAR	0.016238	0.000019
4	DENF-2013-000017	0	2012	51	20/12/2012	28/12/2012	20/12/2012	1/1/2013	534969	33192.61720	36691.75100	2861	TAI SENG	0.016238	0.000019

gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["X output"], df["Y output"]))

gdf=gdf.set_crs(epsg="3414")

gdf.head()

	CASE_ID	CLUSTER_ID	Onset_EYear	Onset_EWeek	Clean_Onset	ONSET_DATE	NEA_ONSET_	NOTIFY_DAT	POSTAL	X output	Y output	cell_id	friendly_n	Shape_Leng	Shape_Area	geometry
0	DENF-2013-000025	0	2012	51	17/12/2012	17/12/2012	17/12/2012	31/12/2012	650536	18857.18765	37512.04877	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18857.188 37512.049)
1	DENF-2013-000018	20130008	2012	51	21/12/2012	31/12/2012	21/12/2012	31/12/2012	536467	32943.97740	36955.30429	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32943.977 36955.304)
2	DENF-2013-000048	0	2012	51	22/12/2012	30/12/2012	22/12/2012	3/1/2013	550421	32663.67780	36774.03736	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32663.678 36774.037)
3	DENF-2013-000109	0	2012	51	22/12/2012	3/1/2013	22/12/2012	4/1/2013	550422	32750.05966	36774.55328	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32750.060 36774.553)
4	DENF-2013-000017	0	2012	51	20/12/2012	28/12/2012	20/12/2012	1/1/2013	534969	33192.61720	36691.75100	2861	TAI SENG	0.016238	0.000019	POINT (33192.617 36691.751)

gdf_latlong = gdf.to_crs(epsg = "4757")

gdf_latlong.head()

	CASE_ID	CLUSTER_ID	Onset_EYear	Onset_EWeek	Clean_Onset	ONSET_DATE	NEA_ONSET_	NOTIFY_DAT	POSTAL	X output	Y output	cell_id	friendly_n	Shape_Leng	Shape_Area	geometry
0	DENF-2013-000025	0	2012	51	17/12/2012	17/12/2012	17/12/2012	31/12/2012	650536	18857.18765	37512.04877	1334	HONG KAH NORTH	0.016238	0.000019	POINT (103.75116 1.35552)
1	DENF-2013-000018	20130008	2012	51	21/12/2012	31/12/2012	21/12/2012	31/12/2012	536467	32943.97740	36955.30429	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (103.87774 1.35048)
2	DENF-2013-000048	0	2012	51	22/12/2012	30/12/2012	22/12/2012	3/1/2013	550421	32663.67780	36774.03736	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (103.87522 1.34885)
3	DENF-2013-000109	0	2012	51	22/12/2012	3/1/2013	22/12/2012	4/1/2013	550422	32750.05966	36774.55328	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (103.87600 1.34885)
4	DENF-2013-000017	0	2012	51	20/12/2012	28/12/2012	20/12/2012	1/1/2013	534969	33192.61720	36691.75100	2861	TAI SENG	0.016238	0.000019	POINT (103.87998 1.34810)

Hexagons Layer

hexagon_path = r"C:\Users\User\Documents\EHI\Projects\VAR Modelling\Data Mgmt\hexagons\hexagons.shp"
hexagons_layer = gpd.read_file(hexagon_path)
hexagons_layer.crs

<Geographic 2D CRS: EPSG:4757>
Name: SVY21
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: Singapore
- bounds: (103.59, 1.13, 104.07, 1.47)
Datum: SVY21
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

hexagons_layer.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
RangeIndex: 4083 entries, 0 to 4082
Data columns (total 11 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    4083 non-null   int64   
 1   TARGET_FID  4083 non-null   int64   
 2   centroid    4083 non-null   object  
 3   cell_id     4083 non-null   object  
 4   friendly_n  4083 non-null   object  
 5   Shape_Leng  4083 non-null   float64 
 6   Shape_Area  4083 non-null   float64 
 7   cell_id_1   4083 non-null   int64   
 8   X           4082 non-null   float64 
 9   Y           4082 non-null   float64 
 10  geometry    4082 non-null   geometry
dtypes: float64(4), geometry(1), int64(3), object(3)
memory usage: 351.0+ KB

hexagons_layer.head()

	OBJECTID	TARGET_FID	centroid	cell_id	friendly_n	Shape_Leng	Shape_Area	cell_id_1	X	Y	geometry
0	1	0	{ "type": "Point", "coordinates": [ 103.674855...	463	GUL BASIN	0.016238	0.000019	463	103.675	0.0	POLYGON ((103.67216 1.29272, 103.67351 1.29507...
1	2	1	{ "type": "Point", "coordinates": [ 103.678899...	514	SAFTI	0.016238	0.000019	514	103.679	0.0	POLYGON ((103.67620 1.33270, 103.67755 1.33505...
2	3	2	{ "type": "Point", "coordinates": [ 103.691040...	634	SHIPYARD	0.016238	0.000019	634	103.691	0.0	POLYGON ((103.68834 1.30213, 103.68969 1.30448...
3	4	3	{ "type": "Point", "coordinates": [ 103.840759...	2410	OXLEY	0.016238	0.000019	2410	103.841	0.0	POLYGON ((103.83806 1.29508, 103.83941 1.29743...
4	5	4	{ "type": "Point", "coordinates": [ 104.018803...	3717	CHANGI BAY	0.016238	0.000019	3717	104.019	0.0	POLYGON ((104.01611 1.31388, 104.01745 1.31623...

hexagons_layer.dropna(inplace=True)

hexagons_layer["cell_id"]=hexagons_layer["cell_id"].astype(int)

hexagons_layer=hexagons_layer[hexagons_layer["cell_id"].isin(df["cell_id"])]

hexagons_layer.info()

<class 'geopandas.geodataframe.GeoDataFrame'>
Int64Index: 1563 entries, 3 to 3742
Data columns (total 11 columns):
 #   Column      Non-Null Count  Dtype   
---  ------      --------------  -----   
 0   OBJECTID    1563 non-null   int64   
 1   TARGET_FID  1563 non-null   int64   
 2   centroid    1563 non-null   object  
 3   cell_id     1563 non-null   int32   
 4   friendly_n  1563 non-null   object  
 5   Shape_Leng  1563 non-null   float64 
 6   Shape_Area  1563 non-null   float64 
 7   cell_id_1   1563 non-null   int64   
 8   X           1563 non-null   float64 
 9   Y           1563 non-null   float64 
 10  geometry    1563 non-null   geometry
dtypes: float64(4), geometry(1), int32(1), int64(3), object(2)
memory usage: 140.4+ KB

hexagons_layer.to_crs(epsg="3414")["geometry"].length

3       1801.332902
20      1801.332975
32      1801.332740
34      1801.332770
76      1801.332915
           ...     
3626    1801.332947
3644    1801.332834
3648    1801.332858
3669    1801.332866
3742    1801.332902
Length: 1563, dtype: float64

hexagons_layer.to_crs(epsg="3414")["geometry"].area/ 10**6

3       0.234173
20      0.234173
32      0.234173
34      0.234173
76      0.234173
          ...   
3626    0.234173
3644    0.234173
3648    0.234173
3669    0.234173
3742    0.234173
Length: 1563, dtype: float64

import folium
from folium.plugins import HeatMap
map_folium = folium.Map([1.35255,103.82580],zoom_start=11.4,control_scale=True)
for _, r in hexagons_layer.iterrows():
    #without simplifying the representation of each borough, the map might not be displayed
    #sim_geo = gpd.GeoSeries(r['geometry'])
    sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
    geo_j = sim_geo.to_json()
    geo_j = folium.GeoJson(data=geo_j,style_function=lambda x: {'fillColor': 'orange'})
   # folium.Popup(r['BoroName']).add_to(geo_j)
    geo_j.add_to(map_folium)
map_folium

Make this Notebook Trusted to load map: File -> Trust Notebook

Plot with 500m grid layer

gdf.head()

	CASE_ID	CLUSTER_ID	Onset_EYear	Onset_EWeek	Clean_Onset	ONSET_DATE	NEA_ONSET_	NOTIFY_DAT	POSTAL	X output	Y output	cell_id	friendly_n	Shape_Leng	Shape_Area	geometry
0	DENF-2013-000025	0	2012	51	17/12/2012	17/12/2012	17/12/2012	31/12/2012	650536	18857.18765	37512.04877	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18857.188 37512.049)
1	DENF-2013-000018	20130008	2012	51	21/12/2012	31/12/2012	21/12/2012	31/12/2012	536467	32943.97740	36955.30429	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32943.977 36955.304)
2	DENF-2013-000048	0	2012	51	22/12/2012	30/12/2012	22/12/2012	3/1/2013	550421	32663.67780	36774.03736	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32663.678 36774.037)
3	DENF-2013-000109	0	2012	51	22/12/2012	3/1/2013	22/12/2012	4/1/2013	550422	32750.05966	36774.55328	2826	UPPER PAYA LEBAR	0.016238	0.000019	POINT (32750.060 36774.553)
4	DENF-2013-000017	0	2012	51	20/12/2012	28/12/2012	20/12/2012	1/1/2013	534969	33192.61720	36691.75100	2861	TAI SENG	0.016238	0.000019	POINT (33192.617 36691.751)

grid_layer.head()

	OBJECTID	Id	Shape_Leng	Shape_Area	geometry
0	1	0	2000.000000	250000.000000	POLYGON ((18162.733 31755.775, 17662.733 31755...
1	2	0	2000.000000	250000.000000	POLYGON ((23162.733 42755.775, 23662.733 42755...
2	3	0	2000.000000	250000.000000	POLYGON ((37662.733 40255.775, 37162.733 40255...
3	4	0	2961.516312	226945.081982	MULTIPOLYGON (((28662.733 48255.775, 28235.903...
4	5	0	1974.913982	248059.096956	POLYGON ((29162.733 28255.775, 29162.733 27755...

gdf_grid = gpd.sjoin(gdf,grid_layer,op="within")

gdf_grid.head()

	CASE_ID	CLUSTER_ID	Onset_EYear	Onset_EWeek	Clean_Onset	ONSET_DATE	NEA_ONSET_	NOTIFY_DAT	POSTAL	X output	...	cell_id	friendly_n	Shape_Leng_left	Shape_Area_left	geometry	index_right	OBJECTID	Id	Shape_Leng_right	Shape_Area_right
0	DENF-2013-000025	0	2012	51	17/12/2012	17/12/2012	17/12/2012	31/12/2012	650536	18857.18765	...	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18857.188 37512.049)	2643	2644	0	2000.0	250000.0
95	DENF-2013-000096	20130036	2013	1	5/1/2013	5/1/2013	5/1/2013	7/1/2013	650533	18705.23952	...	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18705.240 37551.897)	2643	2644	0	2000.0	250000.0
511	DENF-2013-000530	20130036	2013	3	14/1/2013	16/1/2013	14/1/2013	19/1/2013	650532	18695.52780	...	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18695.528 37593.481)	2643	2644	0	2000.0	250000.0
773	DENF-2013-001011	20130087	2013	4	26/1/2013	26/1/2013	26/1/2013	1/2/2013	650535	18772.05585	...	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18772.056 37451.456)	2643	2644	0	2000.0	250000.0
1018	DENF-2013-001029	20130087	2013	5	27/1/2013	27/1/2013	27/1/2013	1/2/2013	650535	18772.05585	...	1334	HONG KAH NORTH	0.016238	0.000019	POINT (18772.056 37451.456)	2643	2644	0	2000.0	250000.0

5 rows × 21 columns

grid_subset = grid_layer[grid_layer["OBJECTID"].isin(gdf_grid["OBJECTID"])]

grid_subset=grid_subset.to_crs(epsg ="4757")

grid_subset.to_crs(epsg = "3414")["geometry"].area/10**6

3       0.226945
4       0.248059
6       0.250000
8       0.250000
13      0.250000
          ...   
3490    0.250000
3491    0.250000
3495    0.250000
3496    0.250000
3497    0.250000
Length: 1482, dtype: float64

import folium
from folium.plugins import HeatMap
map_folium = folium.Map([1.35255,103.82580],zoom_start=11.4,control_scale=True)
for _, r in grid_subset.iterrows():
    #without simplifying the representation of each borough, the map might not be displayed
    #sim_geo = gpd.GeoSeries(r['geometry'])
    sim_geo = gpd.GeoSeries(r['geometry']).simplify(tolerance=0.001)
    geo_j = sim_geo.to_json()
    geo_j = folium.GeoJson(data=geo_j,style_function=lambda x: {'fillColor': 'orange'})
   # folium.Popup(r['BoroName']).add_to(geo_j)
    geo_j.add_to(map_folium)
map_folium

Make this Notebook Trusted to load map: File -> Trust Notebook